#for data upload

rm(list = ls())
setwd(#directory#)
library(readstata13)
library(reshape2)
library(magrittr)
library(dplyr)
library(ggplot2)
data <- read.dta13(#upload_ca.dta#)

View(data)
str(data)

#basic model
f1 <- as.formula(trust_chief_diff~ pro_nsl_protest*factor(pro_establish_after_new)+
                   factor(stance0)+hkecon_diff+welfarefuture_diff+oppompdq_diff)
res1 <- lm(f1,data=data)

summary(res1)[["coefficients"]][, "t value"][c(12,13)] #check the number of coefficient

#create independent variables with different cutoffs
data$pro_nsl_protest_1 <- ifelse(data$nsl_protest1>=1,1,0)
data$pro_nsl_protest_2 <- ifelse(data$nsl_protest1>=2,1,0)
data$pro_nsl_protest_3 <- ifelse(data$nsl_protest1>=3,1,0)
data$pro_nsl_protest_4 <- ifelse(data$nsl_protest1>=4,1,0)
data$pro_nsl_protest_5 <- ifelse(data$nsl_protest1>=5,1,0)
data$pro_nsl_protest_6 <- ifelse(data$nsl_protest1>=6,1,0)
data$pro_nsl_protest_7 <- ifelse(data$nsl_protest1>=7,1,0)
data$pro_nsl_protest_8 <- ifelse(data$nsl_protest1>=8,1,0)
data$pro_nsl_protest_9 <- ifelse(data$nsl_protest1>=9,1,0)

#acgt
data["acgt0"] <- (data["trust_central0"]+data["trust_liberaitonarmy0"]+data["trust_liason0"])/3
data["acgt1"] <- (data["trust_central1"]+data["trust_liberaitonarmy1"]+data["trust_liason1"])/3
data["acgt_diff"] <- data["acgt1"]-data["acgt0"]
a <- "pro_nsl_protest"

ins <- c("chief", "court", "central", "party" ,"legco", "liberaitonarmy" ,"reo" ,"liason","police","acgt")
relab <- function(ins_name){
  if (ins_name=="chief"){
    tl <- "TRUST_CE"
  } else if (ins_name=="court"){
    tl <- "TRUST_COURT"
  } else if (ins_name=="central"){
    tl <- "TRUST_CG"
  } else if (ins_name=="party"){
    tl <- "TRUST_PARTY"
  } else if (ins_name=="legco"){
    tl <- "TRUST_LEGCO"
  } else if (ins_name=="liberaitonarmy"){
    tl <- "TRUST_PLA"
  } else if (ins_name=="reo"){
    tl <- "TRUST_REO"
  } else if (ins_name=="liason"){
    tl <- "TRUST_LO"
  } else if (ins_name=="police"){
    tl <- "TRUST_POL"
  } else if (ins_name=="acgt"){
    tl <- "ACGT"
  }
  return(tl)
}

#different cutoff of nsl_protest1 from 1 to 9
for (i in ins){
  print(i)
  for(k in 3:6){
    b <- paste(a,k,sep = "_")
    if (i != "acgt"){
      dep <- paste("trust",i,"diff", sep= "_")
    } else {
      dep <- paste(i,"diff", sep= "_")
    }
    f <- as.formula(get(dep)~ factor(get(b))*factor(pro_establish_after_new)+
                      factor(stance0)+hkecon_diff+welfarefuture_diff+oppompdq_diff)
    res_k <- lm(f,data=data)
    tvalue <- summary(res_k)[["coefficients"]][, "t value"][c(12,13)]
    if(k==3){
      t_1 <- tvalue[1]
      t_2 <- tvalue[2]
    }else{
      t_1 <- c(t_1,tvalue[1])
      t_2 <- c(t_2,tvalue[2])
    }
  }
  
  df <- data.frame(demo=t_1,esta=t_2,cutoff=c(3:6))
  tlabel <- relab(i)
  
  title_demo <- paste(tlabel,"_DIFF",": CA x PRO_DEMO",sep="")
  file_name_d <- paste(paste("~/Documents/椰林大學/助理/Hong_kong_protest/graph_2cat/robustness/",i,sep=""),"_cutoff_demo.png",sep="")
  png(file=file_name_d, width = 480, height = 360)
  p_demo<-ggplot(data=df, aes(x=cutoff, y=demo, fill=factor(ifelse(cutoff==4, "Model", "Alternatives")))) +
    geom_bar(stat="identity") +
    scale_fill_manual(name = "cutoff", values=c("grey50","red")) +
    labs(title=title_demo, x="Cut Off", y="t-statistics")
  print(p_demo)
  dev.off()
  
  title_esta <- paste(tlabel,"_DIFF",": CA x PRO_ESTAB",sep="")
  file_name_e = paste(paste("~/Documents/椰林大學/助理/Hong_kong_protest/graph_2cat/robustness/",i,sep=""),"_cutoff_estab.png",sep="")
  png(file=file_name_e, width = 480, height = 360)
  p_esta<-ggplot(data=df, aes(x=cutoff, y=esta, fill=factor(ifelse(cutoff==4, "Model", "Alternatives")))) +
    geom_bar(stat="identity") + 
    scale_fill_manual(name = "cutoff", values=c("grey50","red")) +
    labs(title=title_esta, x="Cut Off", y="t-statistics")
  print(p_esta)
  dev.off()
}

#simulation
set.seed(1)
for(j in 1:10){
  t_demo <- rep(NA, 10000)
  t_estab <- rep(NA, 10000)
  data_sim <- data
  
  for(i in 1:10000){
    data_sim$note <- NA
    data_sim <- data_sim[sample(1:nrow(data_sim)),]
    data_sim$note[1:600]<-0
    data_sim$note[601:1225]<-1
    if (j != 10){
      insj <- paste(paste("trust",ins[j],sep = "_"),"diff",sep = "_")
    } else {
      insj <- paste(ins[j],"diff",sep = "_")
    }
    f <- as.formula(get(insj)~ factor(note)*factor(pro_establish_after_new)+
                      factor(stance0)+hkecon_diff+welfarefuture_diff+oppompdq_diff)
    res_k <- lm(f,data=data_sim)
    tvalue <- summary(res_k)[["coefficients"]][, "t value"][c(12,13)]
    t_demo[i] <- tvalue[1]
    t_estab[i] <- tvalue[2]
    if(i %% 1000==0){
      print(i)
    }
  }
  
  #get original value
  f1 <- as.formula(get(insj)~ factor(pro_nsl_protest)*factor(pro_establish_after_new)+
                     factor(stance0)+hkecon_diff+welfarefuture_diff+oppompdq_diff)
  res1 <- lm(f1,data=data)
  t_value <- summary(res1)[["coefficients"]][, "t value"][c(12,13)]
  tbb <- relab(ins[j])
  
  #simulation plot
  demo_graph_file <- paste(paste("~/Documents/椰林大學/助理/Hong_kong_protest/graph_2cat/robustness/",ins[j],sep=""),"_demo.png",sep="")
  png(demo_graph_file, width = 520, height = 360)
  demo_graph_name <- paste("10000 simulations on CA, CA x PRO_DEMO", ": ", tbb, "_DIFF", sep="")
  sup <- t_value[1]
  inf <- -t_value[1]
  if (sup < 0){
    sup <- -sup
    inf <- -inf
  }
  if (ins[j]=="central"||ins[j]=="liberaitonarmy"||ins[j]=="police"||ins[j]=="acgt"){
    hist(t_demo, breaks = 100, xlim=c(inf-0.5,sup+0.5), xlab = "t-statistics", main = demo_graph_name)
    abline(v = t_value[1], col = "red")
  }else{
    hist(t_demo, breaks = 100, xlab = "t-statistics", main = demo_graph_name)
    abline(v = t_value[1], col = "red")
  }
  dev.off()
  estab_graph_file <- paste(paste("~/Documents/椰林大學/助理/Hong_kong_protest/graph_2cat/robustness/",ins[j],sep=""),"_estab.png",sep="")
  png(estab_graph_file, width = 520, height = 360)
  estab_graph_name <- paste("10000 simulations on CA, CA x PRO_ESTAB", ": ", tbb, "_DIFF", sep="")
  hist(t_estab, breaks = 100, xlab = "t-statistics", main = estab_graph_name)
  abline(v = t_value[2], col = "red")
  dev.off()
}




